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The question of the nature of dark matter in the Universe remains one of the 
most outstanding unsolved problems in basic science. One of the best motivated 
particle physics candidates is the lightest supersymmetric particle, assumed to be 
the lightest neutralino. We here describe DarkSUSY, an advanced numerical FOR¬ 
TRAN package for supersymmetric dark matter calculations. With DarkSUSY one 
can: (i) compute masses and compositions of various supersymmetric particles; 
(ii) compute the relic density of the lightest neutralino, using accurate methods 
which include the effects of resonances, pair production thresholds and coannihi¬ 
lations; (iii) check accelerator bounds to identify allowed supersymmetric models; 
and (iv) obtain neutralino detection rates for a variety of detection methods, in¬ 
cluding direct detection and indirect detection through antiprotons, gamma-rays 
and positrons from the Galactic halo or neutrinos from the center of the Earth or 
the Sun. 


1. Introduction 

One of the favourite candidates for the dark matter is a Weakly Interacting 
Massive Particle, a WIMP. In supersymmetric extensions of the standard 
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model, the neutralino emerges as a natural WIMP candidate for the dark 
matter of the universe. 

Over several years, we have developed analytical and numerical tools 
for dealing with the sometimes complex calculations necessary to go from 
given supersymmetric input parameters to actual quantitative predictions 
of the neutralino relic density in the Universe, and of the direct and indirect 
detection rates of neutralino dark matter. In 2000, we made available the 
first public release of DarkSUSY Here we outline the different components 
of DarkSUSY, and comment on the updates and improvements in the coming 
release. The latest version of DarkSUSY can be downloaded from the official 
DarkSUSY website, http://www.physto.se/~edsjo/darksusy/. 

2. Supersymmetric model in DarkSUSY 

DarkSUSY implements the general structure of the Minimal Supersymmet¬ 
ric Standard Model (MSSM), with R-parity and CP conservation (except 
in the quarks CKM matrix). The parameter space of the MSSM is specified 
by 124 a priori free parameters^. 

Most of the DarkSUSY code does not make any specific assumption 
about the 124 parameters (except for R-parity and CP conservation, as 
mentioned). However, in the version of DarkSUSY soon to be released, the 
supersymmetric mass spectrum and the particle mixing matrices are com¬ 
puted under one of two classes of restrictions on the choice of supersymmet¬ 
ric parameters: (1) seven input parameters at the electroweak scale, as was 
available in the previous release (see for details), or (2) five parameters 
at the GUT scale, in the context of minimal supergravity (mSUGRA) im¬ 
plemented via an interface to ISASUGRA^. Beyond these two possibilities, 
the user can in principle implement his own flavor of MSSM. 

DarkSUSY calculates all masses, mixings, and most of the vertices en¬ 
tering the Feynman rules. All of these are available to the user. DarkSUSY 
includes several options for the loop corrections to the Higgs masses®’"^’®’®, 
and for the loop corrections to the neutralino and chargino masses^*^di 

3. Accelerator bounds 

Accelerator bounds can be checked by a call to a subroutine. By modifying 
an option, the user can impose bounds as of different moments in time. 
The default option in the coming public release adopts the 2002 limits by 
the Particle Data Group^^. The user is also free to use his own routine to 
check for experimental bounds, in which case he or she would only need to 
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provide an interface to DarkSUSY. 

4. Calculation of the relic density 

DarkSUSY calculates the neutralino relic density by fully computing the 
thermal average of the effective neutralino annihilation cross section, in¬ 
cluding all resonances, thresholds and coannihilations, and then solves the 
Boltzmann equation numerically with the methods given in 

A major update in the next public release is the inclusion of coanni¬ 
hilations with squarks and sleptons. Coannihilations are important in the 
computation of the relic density when other particles the neutralino can 
convert to are close in mass to the lightest neutralino. These other parti¬ 
cles are abundant at the time of neutralino freeze-out in the early Universe, 
and their annihilation rate must be included in the calculation. The cur¬ 
rent public release of DarkSUSY includes coannihilation processes between 
all neutralinos and charginos; the next release includes coannihilations be¬ 
tween all the above and also squarks and sleptons. 

DarkSUSY includes all two-body hnal states that occur at tree level, and 
neutralino annihilation into gg, 77 and Z 7 that occur at the 1 -loop level. 

5. Detection rates 

The different detection rates for neutralino dark matter have been calcu¬ 
lated by many authors in the past. We will here only give a brief review 
about what is included in DarkSUSY, and which calculations they are based 
on. For a more extensive list of references, we refer to 

5.1. Halo models 

Currently implemented in DarkSUSY is the spherical family of halo profiles 
p{r) oc 1/ ((r/a)'*' [1-1- where e.g. the Navarro, Frenk and 

White prohle^^ is given by (a,/3, 7 ) = (1, 3,1) and the isothermal sphere is 
given by (a,/?, 7 ) = (2,2,0). The velocity distribution is assumed to be a 
standard isotropic gaussian distribution. 

5.2. Direct detection 

These routines calculate the spin-dependent and spin-independent scatter¬ 
ing cross sections on protons and neutrons assuming the quark contributions 
to the nucleon spin from The older set of data from is also available 
as an option, and the user can set his/her own values if desired. 
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5.3. Monte Carlo simulations 

In several of the indirect detection processes described below we need the 
yield of different particles per neutralino annihilation. The hadroniza- 
tion and/or decay of the annihilation products have been simulated with 
Pythia^® and the results have been tabulated. These tables are incorpo¬ 
rated in DarkSUSY. 


5.4. Neutrinos from the Sun and Earth 

Neutralinos can accumulate in the Earth and the Sun where they can an¬ 
nihilate pair-wise producing high energy muon neutrinos. The branching 
ratios for different annihilation channels are calculated and the Pythia 
simulations are used to evaluate the yield of neutrinos. Neutrino interac¬ 
tions in the Sun as well as the charged current neutrino-nucleon interaction 
near the detector are also simulated with Pythia. 

There are routines to calculate a) the neutrino flux, b) the neutrino-to- 
muon conversion rate and c) the neutrino-induced muon flux either differ¬ 
ential in energy and angle or integrated within an angular cone and above 
a given threshold. The new population of neutralinos in the solar system 
described in i9-20,2i from neutralinos that have scattered in the 

outskirts of the Sun) can optionally be included as well. 


5.5. Antiprotons from halo annihilations 

Neutralinos can also annihilate in the Milky Way halo producing e.g. an- 
tiprotons. These propagate in the galaxy before reaching us. We have 
implemented the propagation method described in Optionally, the an¬ 
tiproton fluxes can also be solar modulated with the spherically symmetric 
model of There are also other propagation models^^ available as options. 
The antiproton fluxes are given differential in energy. 


5.6. Positrons from halo annihilations 

Neutralino annihilations in the halo can also produce positrons. The flux of 
positrons is calculated with the propagation model in (with two choices 
of the energy dependence of the diffusion constant). The models in or 
can also be used as an option. The positron fluxes are given differential in 
energy. 
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5.7. Gamma rays from halo annihilations 

Halo annihilations can also produce gamma rays. These are either 
monochromatic (produced from 1 -loop annihilation^® into 77 and .^ 7 ^®) 
or have a continuous energy spectrum (produced from tt® decays in quark 
jets®°). 

The flux of gamma rays can be obtained in any given direction on the sky 
for the user’s choice of (a, (3, 7 ) in the halo profile. There are also routines to 
average the flux over a chosen angular resolution. The continuous gamma 
rays use Pythia simulations to calculate the gamma ray flux (differential 
in energy, or integrated above an energy threshold). 

5.8. Neutrinos from halo annihilations 

Neutralino annihilations in the halo can also produce neutrinos. Although 
the neutrino fluxes are small, DarkSUSY contains routines to calculate (a) 
the neutrino flux, (b) the neutrino-to-muon conversion rate, and (c) the 
neutrino-induced muon flux either differential in energy and angle or inte¬ 
grated within an angular cone and above a given energy threshold. 

6. Conclusions 

Over the years we have developed a numerical package, DarkSUSY, for neu¬ 
tralino dark matter calculations in the Minimal Supersymmetric Standard 
Model, MSSM. Since the initial public release in 2000, the main improve¬ 
ments to the code are the inclusion of mSUGRA models via an interface 
to ISASUGRA and the inclusion of all squark and slepton coannihilations. 
These additions will be available in the next public release. 

We have here presented an overview of what the program can do, and 
refer the reader to the upcoming paper ® and manual where the details 
will be given. A test program, provided with the distribution, shows in 
more detail how DarkSUSY is used. 
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